<html>
	<head>
		<title>Convex Hull in Three Dimensions</title>
	</head>

	<body>
	<script type="text/javascript">
			var result;
			(function () {
			 function generateUnitCirclePoints(n) {
			 var n, a = [], t, p, ct, cp, st, sp;
			 for (i = 0; i < n; i+=1) {
			 t = Math.random()*2*Math.PI;
			 p = Math.random()*Math.PI;
			 ct = Math.cos(t);
			 cp = Math.cos(p);
			 st = Math.sin(t);
			 sp = Math.sin(p);
			 a.push([ct*sp,st*sp,cp]);
			 }
			 return a;
			 }

			 function makeRidge(a, b) {
			 return a < b ? [a, b] : [b, a];
			 }

			 function makeFacet(a, b, c) {
			 return a < b && a < c ? [a, b, c] : b < c ? [b, c, a] : [c, a, b];
			 }

			 function dot3D(a, b) {
				 return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
			 }
			 function vec3D(a, b) {
				 return [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]];
			 }
			 function diff(b, a) {
				 return [b[0] - a[0], b[1] - a[1], b[2] - a[2]];
			 }
			 function add(a, b) {
				 return [a[0] + b[0], a[1] + b[1], a[2] + b[2]];
			 }

			 function normalize(a) {
				 var n = Math.sqrt(dot3D(a,a));
				 if (n === 0) {
					 return a;
				 }
				 return [a[0]/n, a[1]/n, a[2]/n];
			 }

			 function giftWrap(points, a, b, c) {
				 var ut, un;

				 if (!b) {
					 ur = [1,0,0];
					 un = [0,1,0];
				 } else {
					 ur = normalize(diff(b, a));
					 un = normalize(vec3D(c ? diff(c, a) : [1,0,0], ur));
				 }
				 ut = normalize(vec3D(un, ur));

				 function argmin(m, p) {
					 if (p !== a && p !== b && p !== c) {
						 var v = Math.atan2(dot3D(un, diff(p, a)),dot3D(ut, diff(p, a)));
						 return m.min < v ? m : {argmin: p, min: v};
					 } else {
						 return m;
					 }
				 }

				 return points.reduce(argmin, {argmin:[], min:Math.PI}).argmin;

			 }

			 function indexOf(a, n) {
				 var i, j;

				 if (!Array.isArray(a)) {
					 throw new Error('illegal argument');
				 }
				 if (Array.isArray(n)) {
					 for (i = 0; i < a.length; i += 1) {
						 if (!Array.isArray(a[i]) || a[i].length != n.length) {
							 continue;
						 }
						 for (j = 0; j < n.length; j += 1) {
							 if (a[i][j] !== n[j]) {
								 break;
							 }
						 }
						 if (j === n.length) {
							 return i;
						 }
					 }
					 return -1;
				 } 

				 return a.indexOf(n);
			 }

			 function convexHull3D(points) {
				 var a, b, c, cp, cp_v;
				 var ch = { facets : [], ridges : [], vertices : [] };
				 var rq = [];
				 var f, r;

				 function argmin_z(p1, p2) {
					 return p1[2] < p2[2] ? p1 : p2;
				 }

				 a = ch.vertices.push(points.reduce(argmin_z)) - 1;

				 rq.push([ch.vertices.push(giftWrap(points, ch.vertices[a], null, null)) - 1, a]);

				 while (rq.length) {
					 f = rq.pop();

					 r = makeRidge(f[0], f[1]);
					 if (indexOf(ch.ridges, r) >= 0) {
						 continue;
					 }
					 ch.ridges.push(r);

					 cp_v = giftWrap(points, ch.vertices[f[0]], ch.vertices[f[1]], ch.vertices[f[2]]);

					 cp = indexOf(ch.vertices, cp_v);
					 if (cp < 0) {
						 cp = ch.vertices.push(cp_v) - 1;
					 }

					 f = makeFacet(f[0], cp, f[1]);
					 if (indexOf(ch.facets, f) < 0) {
						 ch.facets.push(f);

						 f = f.slice(0)
							 f.unshift(f.pop());
						 rq.push(f);

						 f = f.slice(0)
							 f.unshift(f.pop());
						 rq.push(f);

						 f = f.slice(0)
							 f.unshift(f.pop());
						 rq.push(f);
					 }
				 }

				 return ch;
			 }

			 function facetAngles(ch) {
				 ch.angles = [];
				 var f, a, n, o, p, on, op, facetNormal;

				 for (f = 0; f < ch.facets.length; f += 1) {
					 facetNormal = [0, 0, 0];
					 for (v = 0; v < ch.facets[f].length; v += 1) {
						 n = ch.vertices[ch.facets[f][v]];
						 o = ch.vertices[ch.facets[f][(v+1)% ch.facets[f].length]];
						 p = ch.vertices[ch.facets[f][(v+2)% ch.facets[f].length]];

						 facetNormal = add(facetNormal, vec3D(diff(n, o), diff(p, o)));
					 }
					 facetNormal = normalize(facetNormal);

					 for (v = 0; v < ch.facets[f].length; v += 1) {
						 n = ch.vertices[ch.facets[f][v]];
						 o = ch.vertices[ch.facets[f][(v+1)% ch.facets[f].length]];
						 p = ch.vertices[ch.facets[f][(v+2)% ch.facets[f].length]];

						 ch.angles.push({facet: f, vertex: ch.facets[f][(v+1) % ch.facets[f].length], angle: Math.atan2(dot3D(vec3D(diff(n, o), diff(p, o)), facetNormal), dot3D(diff(n, o), diff(p, o)))});
					 }
				 }

				 ch.angles.sort(function (a, b) {
						 return a.angle > b.angle ? -1 : a.angle < b.angle ? +1 : 0;
						 });

				 return ch;
			 }

			 function cutRidge(ch, r) {
				 console.log(ch, r);

				 ch.ridges = ch.ridges.filter(function (x) { return !(x[0] === r[0] && x[1] === r[1]); });

				 // find face that has only one ridge in ch, and unfold
				 // that, forming a new face, then repeat until no ridges
				 // left, or no such face can be found.
			 }

			 function unfoldOnce(ch) {
				 var f = ch.facets[ch.angles[0].facet];
				 var vi = f.indexOf(ch.angles[0].vertex);

				 cutRidge(ch, makeRidge(f[vi], f[(vi+1) % f.length]));
				 cutRidge(ch, makeRidge(f[vi], f[(vi+f.length-1) % f.length]));
			 }

			 result = facetAngles(convexHull3D(generateUnitCirclePoints(4)));
			 unfoldOnce(result);
			 console.log(result);
			}());
</script>
	</body>
</html>

